// QR.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
//#include <complex.h>//error:不能加h
#include <complex>
#include <math.h>
//#include <E:\projects\QR\QR\Eigen\Dense>
#include <stdio.h>
#include <iostream>
#include <vector>
#include <fstream>
#include <iomanip>
#define N 10
#define E 10e-12
#define L 1000000

//double  Cr[10][10];
//using namespace Eigen;
using namespace std;

//void M_A(Matrix2d &Mk, int r, int m) {
//
//}
double sgn(double x)
{
    double z;
	z = x > 0 ? 1 : -1;
	return z;
}
double fanshu2(double s[10])
{
	double z=0;
	for (int j = 0; j < N; j++)
	{
		z += s[j] * s[j];
	}
	return sqrt(z);
}
void nishangsanjiao(double A[N][N])
{
	setprecision(12); setw(20);//设置输出格式
	//int r = 1;
	double s[10] = { 0 };
	double u[10];
	double h ;
	double c;
	double p[N] = { 0 };
	double q[N] = { 0 };
	double w[N];
	double t = 0;
	//for (int i = 0; i < N; i++)
		for (int r=1; r < N - 1;)
	{
		bool flag = 0;
		for (int l =r+1; l < N; l++)
		{
			if (A[l][r-1] != 0)
			{
				flag = 1;
				break;
			}
		}
		if (flag == 1) 
		{
			for (int j = 0; j < N; j++)
			{
				s[j] = 0;
			}
			for (int j = r; j < 10; j++)
			{
				s[j] = A[j][r-1];
			}
			c = -sgn(s[r])*fanshu2(s);
			//ur
			for (int k = 0; k < N; k++)
			{
				u[k] = s[k] - c*(k == r);
			}
			//pr qr
			h = c*c - c*s[r];
			for (int k = 0; k < N; k++)
			{
				p[k] = { 0 };
				q[k] = { 0 };

			}
			for (int k = 0; k < N; k++)
				for (int l = 0; l < N; l++)
				{
					p[k] += A[l][k] * u[l]/h;
					q[k] += A[k][l] * u[l] / h;
				}
			//t
			t = 0;
			for (int k = 0; k < N; k++)
			{
				t += p[k] * u[k]/h;
			}
			//w
			for (int k = 0; k < N; k++)
			{
				w[k] = q[k] - t*u[k];
			}
			//Ar
			for (int k = 0; k < N; k++)
				for (int l = 0; l < N; l++)
				{
					A[k][l] = A[k][l] - w[k] * u[l] - u[k] * p[l];
				}
			r++;
			if (r == N - 1)
			{
				for (int m = 0; m < N; m++)
				{
					for (int n = 0; n < N; n++)
					{
						//cout << A[m][n] << endl;
						printf("%f  ", A[m][n]);
						//if (n == 4)cout << "\n" << endl;
					}
					//cout << '\n' << endl;
					printf("\n");
				}
			}
			/*
			double H[N][N]=0;
			for (int m = 0; m <= i; m++)
			{
				H[m][m] = 1;
			}
			for (int n = i + 1; n < N; n++)
				for (int m = i + 1; m < N; m++)
				{
					H[n][m]=(m==n)-2*
				}
			*/
		}
		else
		{
			r++;
			if (r == N - 1)
			{
				for (int m = 0; m < N; m++)
				{
					for (int n = 0; n < N; n++)
					{
						cout << A[m][n]<< endl;
						//printf("%f ", A[m][n]);
						if (n == 4)cout << "\n" << endl;
					}
					cout << '\n' << endl;
					//printf("\n");
				}
			}
		}
	
	}

}


void QRfenjie(double A[N][N])
{
	double ur[N], Qr[N][N], wr[N], pr[N];
	for (int i = 0; i < N; i++)
		for (int j = 0; j < N; j++)
		{
			Qr[i][j] = (i == j);
		}

	for (int r = 0; r < N - 1; r++)
	{
		bool flag = 0;
		for (int i = r + 1; i < N; i++)
		{
			if (A[i][r] != 0)
			{
				flag = 1;
				break;
			}
		}
		if (flag == 0)break;
		else
		{
			double dr, cr, hr;
			dr = 0;
			for (int j = r; j < N; j++)
			{
				dr += A[j][r] * A[j][r];
			}
			dr = sqrt(dr);
			//cr
			cr = -sgn(A[r][r])*dr;
			//hr
			hr = cr*cr - cr*A[r][r];

			for (int i = 0; i < N; i++)
			{
				ur[i] = 0;
				wr[i] = 0;
				pr[i] = 0;
			}
			//ur
			for (int j = r; j < N; j++)
			{
				ur[j] = A[j][r] - cr*(r == j);
			}
			//wr,pr
			for (int k = 0; k < N; k++)
			{
				for (int l = 0; l < N; l++)
				{
					wr[k] += Qr[l][k] * ur[l];
					pr[k] += A[l][k] * ur[l];
					//ur[k] += A[k][l] * ur[l];
				}
				wr[k] = wr[k] / hr;
				pr[k] = pr[k] / hr;
				//qr[k] = qr[k] / hr;
			}
			for (int k = 0; k < N; k++)
				for (int l = 0; l < N; l++)
				{
					Qr[k][l] = Qr[k][l] - wr[k] * ur[l];
					A[k][l] = A[k][l] - ur[k] * pr[l];
				}
		}
	}
	/*for (int k = 0; k < N; k++)
	{
		for (int l = 0; l < N; l++)
		{
			cout << Qr[k][l] << " ";
		}
		cout << "\n" << endl;
	}
	for(int k = 0; k < N; k++)
	{
		for (int l = 0; l < N; l++)
		{
			cout << A[k][l] << " ";
		}
		cout << "\n" << endl;
	}*/
	ofstream file_Qr;
	file_Qr.open("output_Qr.txt", ios::out);

	for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < N; j++)
		{
			if (abs(Qr[i][j]) < E)Qr[i][j] = 0;
			file_Qr << setiosflags(ios::scientific) << setprecision(12) << setw(20) << Qr[i][j];// << endl;
																							//cout << A[i][j] << "  " << endl;
																							//		printf("%f ", A[i][j]);
		}
		//cout << '\n' << endl;
		file_Qr<< "\n" << endl;
		//printf("\n");
	}
	ofstream file_Ar;
	file_Ar.open("output_Ar.txt", ios::out);

	for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < N; j++)
		{
			if (abs(A[i][j]) < E)A[i][j] = 0;
			file_Ar << setiosflags(ios::scientific) << setprecision(12) << setw(20) << A[i][j];// << endl;
																							//cout << A[i][j] << "  " << endl;
																							//		printf("%f ", A[i][j]);
		}
		//cout << '\n' << endl;
		file_Ar << "\n" << endl;
		//printf("\n");
	}
//	return 0;
}


int QR(double A[N][N])
{
	//setprecision(12); setw(20);//设置输出格式
	int k = 1;
	int m = N - 1;
	//9中要用的全局
	double qr[N] = { 0 };
	//for (; k < L; )
	for (; m >= 0;)
	{
		//3
		if (m >= 1 && (abs(A[m][m - 1]) < E))
		{
			//A[m][m - 1] = 0;//将小于E的值置零
			printf("one of the eiganvalues is %f\n", A[m][m]);
			m = m - 1;
			//4
			if (m == 0)
			{
				printf("the last of the eiganvalue is %f\n", A[0][0]);
				m = m - 1;
				ofstream file;
				file.open("output.txt", ios::out);

				for (int i = 0; i < N; i++)
				{
					for (int j = 0; j < N; j++)
					{
						if (abs(A[i][j]) < E)A[i][j] = 0;
						file << setiosflags(ios::scientific) << setprecision(12) << setw(20) << A[i][j];// << endl;
						//cout << A[i][j] << "  " << endl;
						//		printf("%f ", A[i][j]);
					}
					//cout << '\n' << endl;
					file << "\n" << endl;
					//printf("\n");
				}
				return 0;
			}
			//else if (m == 0)
			//{
			////	return 0;
			//}
			else
			{
				continue;
			}
		}
		else
		{
			//5
			double s1_real, s1_img, s2_real, s2_img;
			double D[2][2];
			D[0][0] = A[m - 1][m - 1];
			D[0][1] = A[m - 1][m];
			D[1][0] = A[m][m - 1];
			D[1][1] = A[m][m];

			double det = D[0][0] * D[1][1] - D[0][1] * D[1][0];
			double bb = (D[0][0] + D[1][1])*(D[0][0] + D[1][1]);
			double deta = bb - 4 * det;
			double M[N][N] = { 0 };
			//				double I[N][N];
							//I eye Matrix
							/*for (int i = 0; i < N;i++)
								for (int j = 0; j < N; j++)
								{
									I[i][j] = i == j ? 1 : 0;
								}*/
								//5
			if (deta < 0)
			{
				s1_img = sqrt(abs(deta)) / 2;
				s2_img = -s1_img;
				s1_real = (D[0][0] + D[1][1]) / 2;
				s2_real = s1_real;
			}
			else
			{
				s1_real = ((D[0][0] + D[1][1]) + sqrt(deta)) / 2;
				s2_real = ((D[0][0] + D[1][1]) - sqrt(deta)) / 2;
				s1_img = 0;
				s2_img = 0;
			}
			//6
			if (m == 1)
			{
				printf("the last two of the eiganvalues are %f+i*%f,%f+i*%f\n", s1_real, s1_img, s2_real, s2_img);
				std::ofstream file("output.txt");
				file.open("output.txt", ios::out);

				for (int i = 0; i < N; i++)
				{
					for (int j = 0; j < N; j++)
					{
						if (abs(A[i][j]) < E)A[i][j] = 0;
						file << setiosflags(ios::scientific) << setprecision(12) << setw(20) << A[i][j];// << endl;
						//cout << A[i][j] << "  " << endl;
						//printf("%f ", A[i][j]);
					}
					//cout << '\n' << endl;
					file<<"\n"<<endl;
				}
				return 0;
			}
			else if (m >= 2 && (abs(A[m - 1][m - 2]) <= E))//7
			{
			//	A[m - 1][m - 2] = 0;//将小于E的值置零
				printf("two of the eiganvalues are %f+i*%f,%f+i*%f\n", s1_real, s1_img, s2_real, s2_img);
				m = m - 2;
				//4
				if (m == 0)
				{
					printf("the last of the eiganvalue is %f\n", A[0][0]);
					ofstream file;
					file.open("output.txt", ios::out);

					for (int i = 0; i < N; i++)
					{
						for (int j = 0; j < N; j++)
						{
							if (abs(A[i][j]) < E)A[i][j] = 0;
							file << setiosflags(ios::scientific) << setprecision(12) << setw(20) << A[i][j];// << endl;
							//cout << A[i][j] << "  " << endl;
					//		printf("%f ", A[i][j]);
						}
						//cout << '\n' << endl;
						file << "\n" << endl;
						//printf("\n");
					}
					return 0;
				}
				//else if (m == 0)
				//{
				//	return 0;
				//}
				else
				{
					continue;
				}
			}
			else if (k == L)//8
			{
				printf("attaining the most times");
				return 0;
			}
			else
			{
				//9
				//for (int j = 0; j < m; j++)//?????
				double s = A[m - 1][m - 1] + A[m][m];
				double t = A[m - 1][m - 1] * A[m][m] - A[m][m - 1] * A[m - 1][m];
				bool flag_1 = 0;
				double dr = 0;
				double cr = 0;
				double hr = 0;
				double ur[N] = { 0 };
				double vr[N] = { 0 };
				double pr[N] = { 0 }, tr, wr[N] = { 0 };
				double AA[N][N] = { 0 };
				//AA置零
				for (int i = 0; i < N; i++)
					for (int j = 0; j < N; j++)
					{
						AA[i][j] = 0;
					}
				//求Mk
				for (int i = 0; i <= m; i++)
					for (int k = 0; k <= m; k++)
					{
						//A[N][N] = { 0 };
						//for (int j = 0; j < m; j++)
							//for (int l = 0; l < m; l++)
						for (int x = 0; x <= m; x++)
							//for (int y = o; y < m; y++)
						{
							AA[i][k] += A[i][x] * A[x][k];
						}
						M[i][k] = AA[i][k] - s*A[i][k] + t*(i == k);
					}
				//Mk QR分解和Ak+1求解
				//double B[N][N], C[N][N];
				/*for (int i = 0; i < N; i++)
					for (int j = 0; j < N; j++)
					{
						B[i][j] = M[i][j];
						C[i][j] = A[i][j];

					}*/
				for (int r = 0; r <= m - 1; r++)/////////////////////////////////////////////////

				{

					flag_1 = 0;

					for (int j = r + 1; j <= m; j++)
					{
						if (M[j][r] != 0)
						{
							flag_1 = 1;
							break;
						}
					}
					if (flag_1 == 0)
					{
						continue;
					}
					else
					{
						//dr
						dr = 0;
						for (int j = r; j <= m; j++)
						{
							dr += M[j][r] * M[j][r];
						}
						dr = sqrt(dr);
						//cr
						cr = -sgn(M[r][r])*dr;
						//hr
						hr = cr*cr - cr*M[r][r];
						//
						for (int j = 0; j < N; j++)
						{
							ur[j] = 0;
							vr[j] = 0;
							pr[j] = 0;
							qr[j] = 0;
							wr[j] = 0;

						}
						//ur
						for (int j = r; j <= m; j++)
						{
							ur[j] = M[j][r] - cr*(r == j);
						}
						//vr,pr,qr
						for (int k = 0; k <= m; k++)
						{
							for (int l = 0; l <= m; l++)
							{
								vr[k] += M[l][k] * ur[l];
								pr[k] += A[l][k] * ur[l];
								qr[k] += A[k][l] * ur[l];
							}
							vr[k] = vr[k] / hr;
							pr[k] = pr[k] / hr;
							qr[k] = qr[k] / hr;
						}
						//tr

						tr = 0;
						for (int i = 0; i <= m; i++)
						{
							tr += pr[i] * ur[i];
						}
						tr = tr / hr;
						//wr
						for (int i = 0; i <= m; i++)
						{
							wr[i] = qr[i] - tr*ur[i];
						}
						for (int k = 0; k <= m; k++)
							for (int l = 0; l <= m; l++)
							{
								A[k][l] = A[k][l] - wr[k] * ur[l] - ur[k] * pr[l];
							}
						/////////////////////////
						//for (int x = 0; x < m; x++)
						//{
						//	for (int y = 0; y < m; y++)
						//	{
						//		//cout << A[x][y]<<"  " << endl;
						//		printf("%f ", C[x][y]);
						//	}
						//	//cout << '\n' << endl;
						//	printf("\n");
						//}
						////////////////////////////
					}
				}//以上为Mk的QR分解Ak+1求解

			}
			//for (int x = 0; x < m; x++)
			//{
			//	for (int y = 0; y < m; y++)
			//	{
			//		//cout << a[x][y]<<"  " << endl;
			//		printf("%f ", a[x][y]);
			//	}
			//	//cout << '\n' << endl;
			//	printf("\n");
			//}
			k++;
			//	printf("/*********************************%d***************************************/\n", k);
			//	printf("/*********************************%d***************************************/\n", m);

		}

	}
	
}
int main()
{
	
	
	////Matrix<double,10,10> A;
	//for (int i = 0; i < A.rows(); i++)
	//	for (int j = 0; j <A.cols(); j++)
	//	{
	//		if (i != j)
	//			A(i, j)= sin(0.5*(i+1) + 0.2*(j+1));
	//		else 
	//			A(i, j)= 1.52*cos(i+1 + 1.2*(j+1));

	//	}
	/*for (int i = 0; i < A.rows(); i++) {
		cout << A.row(i) << '\n' << endl;
	}*/
	//cout << A.eigenvalues() << endl;
	double A[N][N];
	for (int i = 0; i < N; i++)
		for (int j = 0; j <N; j++)
		{
			if (i != j)
				//A[i][j] = i + j;
			A[i][j]= sin(0.5*(i+1) + 0.2*(j+1));
			else
				A[i][j]= 1.52*cos(i+1 + 1.2*(j+1));
			//A[i][j] = i + 2*j;
		}
	nishangsanjiao(A);
	QRfenjie(A);
	//QR(A);
	getchar();

    return 0;
}